Ab initio energetics and kinetics study of H2 and CH4 in the SI Clathrate Hydrate 
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We present ab initio results at the density functional theory level for the energetics and kinetics 
of H2 and CH4 in the SI clathrate hydrate. Our results complement a recent article by some of the 
authors [G. Roman-Perez et al., Phys. Rev. Lett. 105, 145901 (2010)] in that we show additional 
results of the energy landscape of H2 and CH4 in the various cages of the host material, as well as 
further results for energy barriers for all possible diffusion paths of H2 and CH4 through the water 
framework. We also report structural data of the low-pressure phase SI and the higher-pressure 
phases Sll and SH. 
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Clathrate hydrates are crystalline, ice-like structures 
formed out of water molecules.!^ The water framework 
creates cavities in which gas molecules — typically O2, 
H2, CO2, CH4, Ar, Kr, Xe — can be trapped, which sta- 
bilize the framework. The existence of clathrates was 
first documented in 1810 by Sir Humphry Davy, and 
clathrates became the subject of intensive studies in the 
1930s, when oil companies became aware that clathrates 
can block pipelines.'^ Nowadays, clathrate hydrates are 
of particular interest for two reasons: (i) they are formed 
naturally at the bottom of the ocean, where they are 
often filled with CH4.'^ These deposits mean a tremen- 
dous stock pile of energy, while — at the same time — 
representing a possible global warming catastrophe if re- 
leased uncontrolled into the environment through melt- 
ing; (ii) clathrate hydrates can be used to store H2 in 
its cavities and can be a viable hydrogen-storage ma- 
terial (albeit with moderate hydrogen-storage density).'^ 
For both cases, an understanding of the interaction be- 
tween the guest molecule and the host framework is cru- 
cial for their formation and melting processes, which are 
still understood poorly.!^ In this brief report, we present 
results that elucidate this crucial guest-molecule/host- 
framework interaction and complement a recent paper by 
some of the authors.'^^ We show additional results of the 
energy landscape of H2 and CH4 in the various cages of 
the host material, and we show further results for energy 
barriers for all possible diffusion paths of H2 and CH4 
through the water framework. We also report structural 
data of the phases SI, SII, and SH. 

At low pressure, the methane filled clathrate forms the 
structure SI, consisting of two types of cages. The smaller 
cage is built of water molecules on the vertices of 12 pen- 
tagons with a diameter of approximately 7.86 A,^ and we 
refer to this as 5^^ cage, or alternatively as D cage. The 
larger cage is built of 12 pentagons and two hexagons 
with a diameter of approximately 8.62 A, and we call it 
5^^6^ or T cage. The unitcell has cubic symmetry and 
consists of two 5^^ and six 5^^6^ cages, with a total of 



46 water molecules. At 250 MPa, the structure SI trans- 
forms into a new cubic phase SII, consisting of sixteen 5^^ 
and eight 5^^6^ cages, containing 136 water molecules in 
its unitcell.'^l When the pressure is increased to 600 MPa, 
the structure undergoes another phase transition to the 
hexagonal phase SH.^ This phase has a smaller unit- 
cell of three 5^^, two 4^5^6'^, and one 5^^6^ cages, with 
only 34 water molecules. Very nice graphical representa- 
tions of the different cages and structures can be found 
in Refs. [21 111 El and [8]. While other clathrate- hydrate 
structures exist, structure SI, SII, and SH are the most 
common ones.l^ 

Guest molecules such as H2 and CH4 in the cavities of 
the clathrate hydrates interact with the water framework 
through van der Waals forces. But even the water frame- 
work itself, i.e. the interaction of water molecules through 
hydrogen bonds, has a van der Waals component .l^i To 
capture these effects, we perform here density functional 
theory (DFT) calculations utilizing the truly non-local 
vdW-DF functional, which includes van der Waals inter- 
actions seamlessly into DFT.^iSHIll^g implemented vdW- 
DF using a very efficient FFT formulatiorfl^ into the lat- 
est release of PWSCF, which is a part of the Quantum- 
ESPRESSO package's! For our calculations we used ul- 
trasoft pseudopotentials with a kinetic energy cutoff for 
wave functions and charge densities of 35 Ry and 280 Ry, 
respectively. A self-consistency convergence criterion of 
at least 1 x 10^* Ry was used. All structures were fully 
optimized with respect to volume and atom positions, 
and the force convergence threshold was at least 10"'* 
Ry/a.u. for SI and SH. We have also performed structural 
calculations on SII, but — due to the large unit cell with 
136 water molecules, i.e. 408 atoms — we used a slightly 
less tight force convergence criterium of 5xlO~* Ry/a.u. 
For SI and SH we used a 2 x 2 x 2 Monkhorst-Pack k- 
mesh,!^ while for SII we performed F-point calculations 
only. 

The empty cages are experimentally not stable, but 
they have been shown to be a good starting point for 
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TABLE I. Calculated and experimental lattice constants a and c for the SI, SII, and SH clathrate hydrates. In addition, 
calculated and experimental average nearest-neighbor and next-nearest-neighbor distances are given, as well as bond angles. 
Standard deviations are provided in square brackets. Experimental values for the lattice constants are taken from Ref. [8] 
for methane-filled cages. Experimental values for the averaged quantities are calculated from the structures given in the 
supplemental materials of Ref. [5]. The experimental distances do°_H (rfo^n) seem to be underestimated (overestimated), most 
likely due to the difficulty of accurately determining H positions in X-ray experiments. For SI, neutron scattering experiments 
suggest doLn = 0.97 A and do"-o = 2.755 A.^ Also note that there is some variation in the experimental results for the lattice 
constants in Refs. [Hill and [8]. 

a [A] C [A] d&"_H [A] dgL-H [A] dg,1o [A] Zh-O-H [°] Zo-O-O [°] 

calc. 11.97 — 0.994 [0.001] 1.790 [0.014] 2.781 [0.013] 107.1° [1.0] 108.6° [4.0] 

exp. 11.88 — 0.861 [0.031] 1.911 [0.022] 2.761 [0.017] 109.3° [3.0] 108.7° [3.7] 

calc. 17.35 — 0.994 [0.001] 1.792 [0.016] 2.784 [0.016] 107.1° [0.6] 109.2° [4.3] 

exp. 17.19 — 0.812 [0.016] 1.959 [0.025] 2.768 [0.013] 109.5° [2.0] 109.3° [4.0] 

calc. 12.32 10.01 0.994 [0.001] 1.793 [0.015] 2.782 [0.011] 107.2° [0.9] 108.4° [8.5] 

exp. 12.33 9.92 0.781 [0.040] 1.955 [0.022] 2.775 [0.005] 108.9° [5.1] 108.4° [8.3] 



calculations like ours.l^We have calculated the optimized 
lattice parameters for the SI, SII, and SH structures and 
the results are collected in Table HI We have also cal- 
culated the structures when filled with methane (one 
methane molecule per cage) and filled with hydrogen 
(up to four II2 per cage), but the lattice parameters ex- 
pand less than 0.1% upon filling, such that we have used 
the parameters for the empty cages henceforth. Over- 
all, our optimized lattice constants agree well with pre- 
vious calculation^ and experiment.^ In Table |l] we fur- 
ther analyze the structure of the host materials by cal- 
culating the average nearest-neighbor and next-nearest- 
neighbor distances and important bond angles. In gen- 
eral, the calculated average distances and angles vary 
only insignificantly amongst SI, SII, and SH, whereas 
they show a slightly larger spread for some experimental 
values. As a side-note, for a single water molecule we cal- 
culate do-H = 0.973 A and Zh-o-h = 104.9°, in good 
agreement with the experimental numbers of 0.958 A and 
104. 5° Note that vdW^ is known to give shghtly too 
large binding distances.'^^^ Small deviations are visible 
in the distances do"-H ^^'^ ^o™Hi which in sum mostly 
cancel to give very good agreement with the experimen- 
tal 0-0 distances. Reference [S] also gives the 0-0 dis- 
tances for all structures explicitly as between 2.725 A and 
2.791 A, in remarkable agreement with our calculations. 
Also, our calculated angles Zq-o-o agree very well with 
experiment. However, the good agreement between oxy- 
gen distances and angles — which describe the structure 
as a whole — is closely related to the agreement for the 
lattice constants. 

We next focus on the binding energies of guest 
molecules in the SI structure. In particular, we study 
the binding energies of CH4 and H2 in the D and T 
cages as a function of the number of molecules; results 
are depicted in Fig. [T] for calculations where molecules 
are added to only one cage in the unitcell, while all other 
cages are kept empty. Here, we define the binding energy 
as the energy difference between the "water-framework -|- 
guest-molecules" system minus the energy of the single 
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FIG. 1. Binding energies per molecule for different numbers of 
CII4 and H2 molecules in the D and T cages of SI. While the 
D and T cages can store only one CH4 molecule, the D and T 
cage can store up to two and four H2 molecules, respectively. 



constituents. In case of 71-fold occupied cages, we sub- 
tract n times the energy of the single molecule. Methane 
is a large molecule compared to the cage sizes and it can 
be seen that in both D and T only one methane molecule 
can be stored. Upon adding another methane molecule, 
the binding energy increases drastically. The situation 
is different for the much smaller H2 molecules. In the 
smaller D cage we can store up to two H2 molecules, but 
increasing the number to three or four results in a positive 
binding energy, i.e. work is required to place more than 
two molecules into this cage. Note that the binding en- 
ergy that we find for double H2 occupancy is rather small, 
i.e. —8 meV/per molecule. It is thus likely that at non- 
zero temperatures cages are only singly occupied. Ex- 
perimentally, while the majority of recent work seems to 
favor single occupancy of the D cages (see e.g. Ref. 2Dj), 
there are also reports that propose double occupancy or 
that find inconclusive evidence .l^lHSSl On the other hand, 
the larger T cage can store four H2 molecules. If a fifth 
molecule is added, it escapes through one of the hexago- 
nal faces into the neighboring, empty T cage. Our calcu- 
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FIG. 2. Energy landscape in meV for a rotating methane 
molecule in a D cage. The x and y axes correspond to ro- 
tations about two mutually perpendicular axis. At the (0, 0) 
point of the plot, the hydrogen atoms of the methane point 
exactly toward four oxygens of the D cage. The difference 
between the minimum and maximum energy is 22 meV. 



lated H2 storage capacity of four molecules in the T cages 
is in agreement with experiment .^^ The binding energy for 
one H2 molecule compares well with quantum-chemistry 
calculations on isolated cavities, which give -0.123 eV.I^ 
Overall, our binding energies are slightly smaller than the 
ones in Ref. f6|. Note that quantum motions have been 
neglected in our approach, which may play an important 
role in the binding process and when determining the 
cage occupancy. A more precise treatment requires the 
computation of the corresponding thermodynamic parti- 
tion function, as for example shown in Ref. [26 . Nev- 
ertheless, we consider our calculations for the binding 
energy an important first step that already reveals im- 
portant information. 

It is also interesting to study where and how the H2 and 
CH4 molecules bind in the cages. If only one molecule is 
present in the cages, it binds in the center of the cage. 
Rotations and small displacements of II2 in that situa- 
tion are on an energy scale of approximately 1 meV and 
approach the accuracy of our calculations. At room tem- 
perature, such perturbations are thus easily thermally ac- 
tivated. Since the methane molecule is larger, it cannot 
move/rotate as easily. We have studied the rotation of a 
single methane molecule centered in the D and T cages as 
a function of rotation about two mutually perpendicular 
axes. The energy landscape for this rotation is depicted 
in Fig. [2] for the D cage. The D cage with its twelve 
pentagons and the methane molecule have a related sym- 
metry, which allows us to choose the (0, 0) point of the 
plot such that all methane hydrogens point exactly to an 
oxygen of the host lattice. At this point, hydrogen bonds 
are created and the total energy is the lowest. Upon ro- 
tation of the methane, the hydrogen bonds break and 
the energy increases. The difference between the lowest 
and highest point of this energy landscape is 22 meV, 
suggesting thermal activation of rotations at room tem- 
perature, and quantifying an experimental assumption.'^ 
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FIG. 3. Barriers to diffusion of H2 and CII4 through the water 
framework along different paths, as a function the relative 
progress along the path. The plots are labeled as yl A B, 
where A is the type of the starting cage, B is the type of the 
ending cage, and x refers to either pentagon (p) or hexagon 
(h), indicating the opening being used for traversing. The 
symbols are the results from 12-image NEB calculations and 
the lines are fitted cubic splines, serving as a guide to the eye. 



We have also studied the rotation of a methane molecule 
in a T cage and the results are very similar to the results 
presented in Fig. |2j with the difference that the maximal 
energy barrier is slightly smaller, i.e. 18 meV, not sur- 
prising as the T cage is slightly larger and the methane 
molecule can rotate more easily. Calculations for both 
rotation-energy landscapes have independently also been 
performed using Siesta^'' ^ - and give essentially identical 
results. 

Finally, we present results for barriers to diffusion 
through the water framework in the SI structure. In the 
case of H2 and hydrogen storage this is of much interest, 
as practical storage solutions require fast kinetics, i.e. low 
barriers. In the case of CH4, the barriers can help us un- 
derstand the natural formation of the filled clathrates. 
We have calculated the barriers to diffusion with nudged 
elastic band (NEB) calculations, using 12 images along 
the path from the center of one cage to the center of the 
next cage; the results for all possible paths are plotted 
in Fig. |3] Note that the relaxation of the host lattice 
is crucial to obtain accurate barrier energies,^ and NEB 
calculations allow for such relaxations perpendicular to 
the path automatically. The plots are labeled as A A B, 
where A is the type of the starting cage, B is the type 
of the ending cage, and x refers to either pentagon (p) 
or hexagon (h), indicating the opening being used for 
traversing. Note that for the path T A D there is only 
one choice of opening, i.e. a pentagon, and the path is 
not symmetric as the distance from the center of T to 
its edge is longer than the corresponding distance in the 
D cage. Furthermore, this path's end energy is differ- 
ent from its starting energy, since the guest molecules 
are binding with different binding energies in the cages 
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T and D, as already evident from Fig. [T] The lowest 
barriers for H2 and C H4 diffusion agree well with previ- 
ous calculations But, our H2 diffusion barrier is an 
overestimation with respect to a recent NMR experiment, 
which gives 0.03 eV and warrants further investigation.^- 
The barriers are in general smaller for diffusion through 
hexagons, simply because these openings are larger. 

For hydrogen-storage applications, the low barrier of 
~0.3 eV between T cages (going through a hexagon) is 
important. Through these T-cage channels, which thread 
through the material in all three dimensions, the hydro- 
gen can quickly be absorbed or released. However, to 
achieve the material's full storage potential, some hy- 
drogen molecules will also have to get into the D cages, 
with a much higher barrier of ~0.75 eV. The large bar- 
rier of ^1.4 eV for methane diffusion suggests that the 
methane molecules get trapped while the clathrate is 
formed, rather than diffusing into an already existing 
empty clathrate. 



To conclude, we have performed an ab initio study of 
structural, energetic, and kinetic properties of the guest 
molecules H2 and CH4 in hydrate clathrates. We have 
also shown first results for the difRcult-to-modcl, high- 
pressure phase SII with a large unit cell, finding good 
agreement with experiment. While we have used vdW- 
DF for our study, it is conceivable that its successor, 
vdW-DF2,^^ may further improve upon our results. We 
encourage additional studies of the hydrate clathrates us- 
ing vdW-DF2, also including other types of cages, and 
more detailed studies of the SII phase, which is one of the 
more promising phases amongst the hydrate clathrates 
for hydrogen-storage applications. 

We would like to dedicate this report to the memory of 
Prof. David Langreth, who passed away just days before 
it was submitted — he is the "father" of vdW-DF and his 
research inspired many. All calculations were performed 
on the WFU DEAC cluster. 
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